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(N 1 ABSTRACT 



Context. The strong electric fields associated with magnetic reconnection in solar flares are a plausible mechanism to 
accelerate populations of high energy, non-thermal particles. One such reconnection scenario, in a fully 3D geometry, 
occurs at a magnetic null point. Here, global plasma motion can give rise to strong currents in the spine axis or fan 
plane. 

Aims. To understand the mechanism of charged particle energy gain in both the external drift region and the diffusion 
region associated with 3D magnetic reconnection. In doing so we aim to evaluate the efficiency of resistive spine and 
fan models for particle acceleration, and find possible observables for each. 

Methods. We use a full orbit test particle approach to study proton trajectories within electromagnetic fields that are 
exact solutions to the steady and incompressible magnetohydrodynamic equations. We study the acceleration physics 
of single particle trajectories and find energy spectra from many particle simulations. The scaling properties of the 
accelerated particles with respect to field and plasma parameters is investigated. 
Results. For fan reconnection, strong non-uniform electric drift streamlines can accelerate the bulk of the test particles. 
The highest energy gain is for particles that enter the current sheet, where an increasing "guide field" stabilises particles 
against ejection. The energy is only limited by the total electric potential energy difference across the fan current sheet. 
The spine model has both slow external electric drift speed and weak energy gain for particles reaching the current 
O ' sheet. 

Conclusions. The electromagnetic fields of fan reconnection can accelerate protons to the high energies observed in solar 
flares, gaining up to 0.1 GeV for anomalous values of resistivity. However, the spine model, which gave a harder energy 
spectrum in the ideal case, is not an efficient accelerator after pressure constraints in the resistive model are included. 
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1. Introduction electrons is a significant fraction of the emission site den- 
sity, setting tough efficiency constraints on any proposed 

Observations of Hard X-ray (HXR) and y-r ay emission from acceleration mechanism. 

solar flares by the RHESSI space telescope (|Lin et al-llgfJOj) ft i s we ll accepted that magnetic reconnection plays 

^ ; suggest that a large proportion of magnetic energy is con- the key role in the dissipation of magnetic energy during 

" verted into kinetic energy of non-thermal accelerated par- a g are anc j t h ere is a growi ng body of obse rvational sig- 

ticles. The dominant HXR sources are chromospheric foot- na t U res for the process fsee lMcKenzi3l201lh . The site of 

points of the flaring loops, at which there is continuum reconnection in the standard (CSHKP) flare model (eg. 

KJ . free-free emission from a be am of energet ic electrons in col- iPriesti I2000T) is above the thermal loops, not in disagree- 

•jjj ■ lision with ambient plasma (jBrown|[197l|) . This continuum ment wit i 1 tne site of the non -thermal coronal HXR source. 

03 \ spectrum gives beam electro n ener gies from around 10 keV Super-Dreicer (|Dreicerlll959l) electric fields associated with 

up to almost 100 MeV dLj^|2006(). There is line emission reconnection are one plausible mechanism for particle ac- 

at the y-ray end of the spectrum from processes involv- ce leration and much theoretical work has been done to in- 

ing accelerated ions such as neutron-ca pture and nuclear vestigate the efficien cy of this mechanism (for review, see 

de-excitation (see eg. IVilmer et a Ll l2011l for a review), sug- IZharkova et~aHl201lh . 

gcsting ions with energies up to ~ 100 MeV/nucleon. ^ , , , , . , „ . , . 

00 of i Early work on charged particle trajectories within a re- 

When the emission from the foot-points is weak, or connecting current sheet concentrated on single particle 
when they are occulted by the solar limb, a weaker HXR motion and e nergy gain in idealised field configurations, 
emission so urce is sometimes o bserved above the top of the ISpeiserl (|l965t ) found that charged particles in the simplest 
flare loops (|Masuda et al.fl994l) . Recent observations of two current sheet, of oppositely directed magnetic field and con- 
such flares indicate that this emission is non-thermal and stant electric field, are trapped so the energy gain is limited 
that the source is actually the accele ration site for a sig- only by the sheet length. With an additional small and con- 
nificant number of en ergetic electrons (|Krucker et al.ll2~010t stant magnetic field component normal to the current sheet 
llshikawa et al.ll201ll ). The estimated number of energetic plane, the particles are turned by 90° and ejected from the 
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current sheet into the external drift region. Zhu k Parks 
(|l993t ) JLitvinenko k Somovl (|l993h and ILitvinenkol (jl996ft 
also included a third component of the magnetic field, a 
guide field parallel to the electric field. Above a critical 
guide field the trajectory is stabilised against ejection and 
the ene rgy gain is once again only bounded by the sheet 
length (|Litvinenkolll996tl . 

This early analytical work was extended with 2D (or 
2.5D where the fields are invariant in the third dimen- 
sion) test particle simulations, many of which have been 
carried out using simple prescribed magnetic and electric 
fields that would be expected in a reconnection solution 
to the MHD equations. These simulations consider the ef- 
fect of the guide field on traject ories and energy spectra 
withi n the Harris current sheet (IZharkova k Gordovskvvl 
12004 120051: IWood k Neukirchl 120051) and magnetic X 
points (IVekstein k 3rownind Il997t lHannah k Fletcher! 
l2006yHamilton et al.ll200jft. 



J eerikhuisen et al. (12002ft and ICraig fc ljitvincnko 
|) used magneti c and electric fields from the exact an- 
alytical solutions of ICraig k Hentonl (|1995l ) to the 2D in- 
compressible, resistive MHD equations. Also, an approach 
combining numerical MHD simulations with a test parti- 
cle co de has also been used to study 2D forced reconnec- 
tion (jGordovskvv et all l2010alfbT ). These simulations can 
include compressibility and time evolution, making them 
more realistic for coronal plasma. However, analytical solu- 
tions are essential to study acceleration due to reconnection 
in very large Lundquist number plasma at present . 

The complexity of the coronal magnetic field in a flar- 
ing Active Region motivates the study of test particle 
motion in fully 3D reconnection geometries. Reconnection 
models in 3D are comparatively new, but it is clear that 
there are significant qu alitative diffe rences from the famil- 
iar 2D models (see eg. iPontinlfeOllft . Reconnection in 3D 
can occur both with and without magnetic null points. 
However, the simplest 3D magnetic configuration is based 
on the potential magnetic field about such a null point. 
Here, the solenoidal condition defines the magnetic topol- 
ogy of a ID spine lin e and a 2D fan su rface (called 7- 
line and E-surface by lLau k Finn! Il990ft that separates 
different magnetic flux domains (a linear description of 
magn etic configurations at nulls is given by iParnell et al.l 
1996). Although these null points cannot be measured in 
the corona at present there is some indirect evidence for 
their existence. Nulls are common features of magnetic 
topology models that reconstruct the magnetic field from 
photo spheric magnetograms into the corona (see Longcopc 
2005, for a review). The application of this method at 
several flare sites suggests the importanc e of these nulls 
in certain flares dDes Jardins et al.l l2009t lAulanier et al.l 
l2000t iFletcher et all 2001ft . Reconnection at magnetic null 
points is also thought to be important for the emer- 
gence o f new flux from beneath the photosphere into the 
coron a (|T5rok et al.l 120091: iMaclean et all l2009t iLiu et al.l 

Urn!. 

The type of reconnection that occurs at a 3D null de- 
pend s upon the magnetic c onfiguration and global plasma 
flow. iPriest k Titovl (JT996) proposed two models of recon- 
nection using a potential magnetic field and prescribed 
global flows that satisfy the ideal MHD equations. In ideal 
spine reconnection a shear flow across the fan plane causes 
frozen-in flux inflow that converges on the spine axis. At 
the spine the field reconnects in the presence of singular 



electric field. For ideal fan reconnection the singular elec- 
tric field occurs in th e fan plane, driven b y a shearing flow 
across the spine axis. ICraig et all (|1995l) . ICraig k Fablind 
(|l996ft and lCraig et al.1 (|l997ft found exact solutions to the 
steady and incompressible resistive MHD equations at 3D 
null points by considering a flux-pileup disturbance field 
superposed with the background potential magnetic field. 
The disturbance field induces a current sheet in the spine 
axis and in the fan plane for the resistive spine a nd resistive 
fan models respectively. ICraig k Fablind (|l998ft found cor- 
responding time-dependent solution s to these steady mod- 
els, an d the numerical simulations of iHeerikhuisen k Craid 
(|2004ft found reconnection scalings in agreement with both 
steady state and time dependent models at peak recon- 
nection rate. The 3D MHD simulations of iPontin et al.l 
(|2007al fbh also found a hybrid of the spine and fan mod- 
els, named spine-fan reconnection, when compressibility is 
included. Recent numerical and analytical study gives addi- 
tional models for null reconnection when the global plasma 
motion is rotational rather than a shear fl ow (see for review 
IPriest k Pontinll2009t IPontin et al.ll201lD . 

It is not yet known if these nul l points are effective parti- 
cle accelera tors. Previous work bylDalla k Browning] (j2005, 
I2006LI2008D and lBrowning et alJ (l2010h use d theldeal elec- 
tromagnetic fields of Priest k Titovl (|l996ft in a test parti- 
cle code, finding the ideal spine reconnection model was 
effective to accelerate protons and electrons to high en- 
ergies (max ~ 10 7 eV) for solar coronal parameters. The 
ideal fan reconnection model was less effective for protons, 
partly as the geometry of the electric drift streamlines was 
less efficie nt at deliv e ring p articles to regions of high elec- 
tric field. iGuo et all J2010;) used null point magnetic and 
electric field configurations from MHD simulations, finding 
that strong electric fields due to convective plasma mo- 
tion can be efficient at accele rating protons but less so 
for electrons. Litvinenkol (|2006ft used the WKB method of 
iBulanov k Cad (11988ft to show that single protons and elec- 
trons close to the null in the reconnecting fan current sheet 
can achieve the high energies observed in flares. However, 
this energy is limited as the particles become unstable in 
the sheet, due to the potential background field, and are 
ejected. 

In this paper we examine test particle trajectories and 
energy spectra of protons in electromagnetic fields which 
are exact solutions to the 3D, steady-state, incompress- 
ible and resistive MHD equations at magnetic null points. 
These ar e the resist i ve spi n e and resistive fan solu tions 
given in ICraig et al.l ((l995h : ICraig k Fabling! {1996) and 
ICraig et al.1 ()1997ft . We consider trajectories that start both 
in the outer ideal region, for compar ison with particle acce l- 
eration results in the ideal models (jBrowning et al.l [20101. 
and those that start directly inside the resistive fan current 
sheet, to compare with the analytical work of ILitvinenkol 

(pool . 

The paper is organised as follows. In Section [5] we de- 
scribe the model fields used, along with parameters cho- 
se n considerin g the p ress ure constraints and opti misations 
of ICraig et all (|1997f) and lCraig k Watsonl {2000). We also 
derive the electric fields and potentials used in the code, 
and give approximate external drift velocity scalings. In 
Section [3] we give the results of test particle simulations for 
thermal distributions of protons starting in the drift region 
for each model. We choose typical high energy particles 
from these simulations and follow the single trajectories to 
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understand the energy gain mechanism. We see how drift 
times and energy spectra scale with the different parameter 
choices in the fan model. In Section U we give a summary 
and conclusions on the efficiency of each model for acceler- 
ating protons. 

2. Model Fields and the Test Particle Code 

We consider two of the reconne c tion solutions at 3D 
nu ll points found b y ICraig et al.l (119951) and dev eloped 
in ICraig fc Fablmel (11996ft and ICraig et al.l ffl997h . The 
fields satisfy the resistive, steady-state, and incompressible 
MHD equations. These are normalised in the usual way by 
the characteristic magnetic field strength Bq, at a typical 
length scale Lq from the null point, and by density pq. This 
choice leads to the natural units for velocities in terms of 
the external Alfven speed, vq = va& — Ba/ \ZPolJ-o- The 
thermal pressure is normalised by the dynamic pressure at 
Alfven speed, po = Pav 2 Ae , so that the dimensionless pres- 
sure on the Lq boundary is half the plasma beta p e = /3 e /2. 
The dimensionless resistivity, 77, is given by 



V 



L VAe^0 



S~ 



(1) 



where r\ d is the dimensional resistivity (Spitzer resistivity 
in the case of purely collisional plasma), po is the mag- 
netic permeability, and S is the Lunquist number which is 
typically very large in the solar corona S ~ 10 12 — 10 14 . 

For completen ess, the main propertie s of the solution 
are giv en here; see ICraig k. Fabl ing (199<|) and ICraig et al.l 
(1997)) for more detail. After normalisation the govern- 
ing equations consist of the momentum equation, which in 
curled form is 

(u • V)w - (w • V)u = (B ■ V) J - ( J • V)B, (2) 

and the induction equation, 

(u ■ V)B - (B ■ V)u = r)V 2 B, (3) 

with the solenoidal and incompressibility conditions, 

V-B=0, V-u = 0. (4) 

Here J is the current density and u> is the vorticity in 
terms of the bulk plasma velocity u. In this normalised form 
they are 

J = V x B, uj = V x u. (5) 

The three dimensional analytic solutions of ICraia: et al.l 
(fl99lLlc7aTg fc Fabling|(fT99l and lcTaig et al.l (119971) have 
magnetic and flow fields of the form 

B = \P + Q, (6) 

u = P + XQ, (7) 

where the scalar < A < 1 gives the shear between the 
B and u fields. The vector field P(x,y,z) is a potential 
background field of strength a, and Q is a disturbance field 
of strength B s which gives rise to current in the models. 

For comparison wit h the particle acceleration results at 
3D nulls in ideal MHD (|Dalla fc Browning|l2005l) we choose 
the z-axis to be aligned with the spine, with z = as the 
fan plane. It must be noted that this choice of axis differs 
from that used bv ICraig et al.l (|1997[ ). We study only the 



proper radial null (jPriest fc Titovlll996|) where the back- 
ground magnetic field lines in the fan plane lie in the radial 
direction. This background field is then written as 



P = -(xx + yy -2zz), 



(8) 



with a giving the sign and strength of the field. For the 
spine model, the displacement field distorts the fan plane 
in the z-direction Q s — Z(x,y)z. For the fan model, it 
distorts the spine in the x-d irection Q P = X (z)x (the 
more general fan case given in ICraig et al.l (|1997l) of Q F = 
X(z)x + Y(z)y has not been covered here). 

2.1. Spine Analytic Fields 

The disturbance field for the spine model in cylindrical co- 
ordinates (r, </>, z) is 

Q s - Z(r, cj>)z = -£—M ( |, 2, ~) sin(0)£, (9) 




(|Craig et al.l fl997) in terms of the confluent hypergeomet- 
ric (K ummer) function M(a,b,() (|Abramowitz fc Stegunl 
1 19T2T ) . The flux pile- up factor B s gives the approximate 
strength of the magnetic field at a dimensionless distance 
r n from the spine axis, where r v is defined as 



(10) 



It is the radius of a cylindrical region centred on the spine 
axis w here resistive effects become significant (jCraig et al.l 
fl997h . 

The form of the displacement field in equation © is only 
a solution to the governing equations provided a < 0. This 
gives frozen-in plasma inflow along the fan plane converging 
on the spine and outflow in the ±z directions away from 
the null point. The magnetic field in the outer (ideal) region 
is also directed inwards along the fan plane and outwards 
along the spine axis. S ome representative magnetic field 
lines are shown in Figure [Jau the displacement term shears 
the fan plane at <j> = ±tt/ 2 towards the spine axis while the 
fieldlines in <f> = 0, it of the fan plane remain perpendicular 
to the spine. 

To integrate particle trajectories using a test particle 
model we require the electric field. We calculate this from 
the uncurled form of equation ([3]) as 



E(r, 



r\dZ „ 
r 

r dcf> 



<P, (11) 



where P r — ar/2 is the radial part of the potential field. 

This electric field is curl-free (as required for steady- 
state) and so we can calculate the electric potential V to 
use as a check of energy conservation. This can be found 
by integrating E = -VF to get 



V(r, (f>) = cos</> 



a (I- A 2 )r 2 /(r) 



rjrf(r) 



(12) 



where f(r) is the radial part of the displacement field in 
Z{r,(j>) = f(r)$m(j>. 

We can study the behaviour of these fields at small 
and large distances using the truncated power series and 
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asymptotic formulae for the K ummer function respectively 
(| Abramowitz fe Steeunl fl972l) . For all our cases the third 
argument in the Kummer function is negative. We have for 
< £ « 1, 

M(o,6,-0w (13) 

and for £ 3> 1, 



M(a,b,-0 



T(b) 



T(b- 



(14) 



in terms of the Gamma function T(b). Near the spine axis 
r f=a 0, the only contribution to the electric field is from 
current in the x-direction and 



E(r « r„) w ??J(0) 



(15) 




The full current distribution is plotted in Figure [ J^b)| it 
forms two cylindrical vortex structures that are localised 
with respect to the resistive region and invariant in the z- 
direction. 

At large distances from the spine, the electric field goes 

as 

—1r\B s sin t 



(a) 



E(r » r„ 



-0. 



(16) 



This h as the same functiona l form as the ideal spine solu- 
tion of lPriest fc Titovl (|l996h . the subje ct of previous work 
on 3D null-p oint particle accelera t ion bylDalla fe B rowning 
(|2005l I2006h and Browning et all (|2010l) . Indeed, by simple 



choice of parameters we could set the magnetic and electric 
fields to asymptotically ma tch those of the ideal case . Some 
care must be taken here as iDalla fc Browning] (|2005| ) stud- 
ied positive nulls, where a > 0, and with E(0 < <j> < w) > 0. 
We have opposite sign for both electric and magnetic fields 
giving the same electric drift inflow quadrants but different 
sign for the convective electric field. Particles that become 
non-adiabatic in the external region r S> r n , and gain en- 
ergy parallel the electric field, wil l rotate about the spine in 
the o pposite direction to those in lDalla fc Browning! (|2005l 
[2006h . In this paper we will only qualitatively compare par- 
ticle trajectories in the ideal and resistive spine models as 
an asymptotic match will give rise to unphysical hydromag- 
netic pressures on the edge of the resistive region r ~ r v 
that were absent in the simplified ideal model (see below). 

The thermal pressure profile for the spine model can be 
found from i ntegrating the unc urled form of equation ©. 
It is given in ICraig et all (|1997| ) to be 



2 
1 

-1 

-2 



x/r 



(b) 

Fig. 1: a) Representative magnetic field lines for the spine 
model with parameters A = 0.75, B s = 3.4, a — — 2, r\ — 3 x 
1CP 3 . The field lines are seeded from the top and base of the 
spine axis, b) Showing the direction and relative strength of 
the current in a plane of constant z for the same parameters. 
Here, r v = v / 4?y sa 0.12 is the size of the resistive region 
centred on the spine axis. 



P = P0 --{P 2 + Z 2 ) 



XazZ, 



(17) 



where po is the gas pressure at the null point, the first term 
inside the brackets is due to dynamic pressure from the 
background flow and the other two terms are from balance 
with magnetic pressure. All terms except for po are nega- 
tive, as a < 0, so constraints must be put on the values of a 
and B s in or der to avoid unphysical n egative pressures as 
discus s ed in iLitvinenko et all (119961); iLitvinenko fc CraigJ 
(|1999h : ICraig et all (|1997T ) and lCraig fc Watsonl (|2000H . We 
give some of the arguments here for the sake of complete- 
ness (see above references for more detail). 

The strong electric field (fast electri c drift) simulations 
for th e ideal spine model studied by IDalla fc Browning] 
(f2005h were characterised by the dimensional value of the 
electric field E = 1500 V/m on the r = 1, (f> = n/2 



boundary (or normalising by suitable solar coronal values, 
v Ae = 6.5 x 10 6 ms- 1 , B Q = 0.01 T, gives E « 1/40). 
This value can be equated with equation (JTBJ) . Crucially, 
to match the external electric field in the resistive model 
to the fixed amplitude electric field in the ideal spine re- 
connection model requires the scaling B s ~ r?" 1 as 77 is re- 
duced to suitable solar coronal values (jCraig et all (|l997l) 
showed that if we require displacement field at the bound- 
ary Z(l,7r/2) ~ 1 this also gives B s ~ i]^ 1 ). However, this 
scaling gives rise to large magnetic pressure on the sheet 
edge. The maximum of the displacement field occurs at 
where Z(r n ) ~ B s giving magnetic pressure there 



from equation ((HJ) Z 2 - B 2 S 



n 



To avoid negative 



thermal pressure in the model this requires the null point 
pressure po > (Z(r v )) 2 ~ rf~ 2 which is unphysically large 
for the values of r\ considered. 
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ICraig et al.l JT997) showed that B s must be limited to 
a saturation value on r = r n , giving weak electric fields 
and small amplitude displacement field on the boundary 
Z(l) <C 1. Also, at r — 1, z <C 1 we have dynamic pressure 
due to bulk fluid inflow p ps p — P 2 where P(l) ~ a. 
We must constrain a < B s or this dynamic pressure will 
require the gas pressure at the null to be even larger. The 
maximum value we can take for po is the largest possible 
external hydromagnetic pres sure p e .max availab le to drive 
the reconnection. We follow ICraig et al.l (|1997( ) and take 
= B emax/ 2 where the maximum external magnetic 
field is that of a sunspot at the photosphere, B e ^ max = 0.3 
T. This gives a normalised saturation value B s ^ max = 30. 

So far we do not know the value a should take, but 
expect that the bulk fluid exhaust from the reconnection 
region is of the order of the local Alfven speed. The exhaust 
on the edge of the current sheet at a global distance from 
the null, r = r v , 4> = 7r/2, z = 1, is given by 



|«(r„,f,l)|«Afl a 



where the local Alfven speed is 

\v A {r v ,%,l)\ = \B{r n ,z,l)\^ B s - \a 

for our choice of normalisation. As we are not interested 
in the case where A = 1 (where there is no shear between 
the velocity and magnetic fields) we have a ~ — B s for 
Alfvcnic exhaust. This is the largest magnitude of a we 
can take without having problems due to dynamic pressure. 
It also leads to the thinnest current sheet and thus max- 
i mises the current density in the resistive region. However, 
as Craig & Watson (2000) show, the dissipation rate is 



W V =1] / J 2 dV « irr)B s 



(18) 



which has no a dependence as the increase in current den- 
sity due to resistive region thinning is cancelled by the r 2 
dependence of the total dissipation volume. 

The electric drift velocity in the external region is given 

by 



v E (r > r v ) 



-2zf 



r]B s sin c 
\\a\y/ir \r(r 2 /A + z 2 ) 



[VAe 



(19) 



which is very slow when \a\ = B s . It is thus necessary to 
limit the magnitude of a so that results can be obtained 
with reasonable integration times. For the simulations in 
Section [3] we use B s = 10, a = —0.1, this limits the re- 
connection exhaust close to the spine current sheet to sub- 
Alfvenic speeds. 

2.2. Fan Analytic Fields 

The displacement field for the fan model is 



Q F = X(z)x 



B,z /3 3 -z 2 
-M 



V 



i/2- 



4' 2' 2fj 



x. 



(jCraig et al.lll997l ). We define z v as 

z n = \, 



2<n 



a(l - A 2 ) ' 



(20) 



(21) 



the approximate height at which X takes the maximum 
value, X max £3 B s . It is a measure of the height of a resistive 



z/L n 




Fig. 2: Representative magnetic field lines for the fan solu- 
tion with parameters A = 0.75, B s = a = 10, rj = 10 -6 . 
The lines are again seeded from the top (solid lines) and 
base (dashed lines) of the spine. 

region centred on the fan plane, z = 0. This form of solution 
is only valid for a > which gives a positive null point, 
the field is washed in from the global boundaries at z = 
±1 and it exits the simulation box radially along the fan 
plane. Some representative magnetic field lines are shown 
in Figure the displacement field shears the spine axis as 
it approaches the fan plane giving rise to strong current 
inside the resistive region. 
The electric field is 

E = y [ V X'(z) - (1 - X 2 )P z X(z)] +z[(l- \ 2 )PyX{z)] , 

(22) 

and the electric potential is 

V(y, z) = -ay(l - A 2 ) [fjX'(z) + zX{z)] , (23) 

where J(z = 0) = X'(0)y is current density at the centre of 
the sheet. The current only has z-dependence; it is infinite 
in extent in the x and y directions. This is cle arly unrealis- 
tic, alt hough resistive MHD simulations by iPontin et al.l 
(j2007bl) find that spine-fan reconnecting current sheets 
formed due to shear flows around a null point spread out 
along the fan plane in the incompressi ble limit. Note that 
analytic multiple null solutions found bv lCraig et al. (1999) 
have finite current sheets, avoiding this problem. In our sim- 
ulations below we consider particle acceleration only within 
a restricted range of 5 Lq, effectively limiting the size of the 
current sheet. 

The thermal p ressure profile for the fan model 
is (jCraig et al.lll997[ ). 



p = p Q -(P 2 + X 2 )/2 - aXxX/2. 



(24) 



However, in this case a displacement field of order unity 
on the z = 1 boundary, X(l) ~ 1, gives the scaling 
B s ~ T? -1 / 4 (jCraig et al.lll997|) . This gives much weaker hy- 
dromagnetic pressure on the current sheet edge compared 
to the spine model but it is still too large for the values of 
rj considered. Again we saturate B s ^ max = 30 and we have 
a < B s to avoid problems from dynamic pressure. 

I Craig fc Watsonl (|2000h show that the Ohmic dissipa- 
tion rate from the fan model is 



W n = T) J J 2 dV 



riB 2 / Zri 



(25) 
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and so in this case, for fixed (saturated) B s , the maximum 
dissipation occurs with the thinnest current sheet (the so 
called optimised solution). The thinnest sheet we can have 
subject to the dynamic pressure constraint is when a = B s 
(for any fixed value of A). Also, as this choice gives the 
largest current density, it maximises the resistive electric 
field within the sheet which is interesting for particle ac- 
celeration. As above, this choice of a sets the bulk fluid 
exhaust at x 2 + y 2 = 1, z = z v to the local Alfven speed. 

Using the asymptotic approximation (|14p we find that 
the z-component of the electric drift that brings the parti- 
cles to the fan plane is, for x, z ^> 77 1 / 2 , 



VEz 



(1 - x 2 )p x p z x{z) ^ b3/4 



XP 2 



1/4 



(26) 



for the optimised solution a — B s . This gives electric drift 
inflow for positive x, z (as P z < 0), and outflow for positive 
z and negative x. It is much faster than the spine case due 
to the more favourable scaling with resistivity. There are 
also fast drift streamlines in the x-y plane that that can 
be found from the numerical (or approximate analytical) 
solution of 

dx dy ,„„, 
= — , (27) 

VEX V E y 

we numerically plot these streamlines on top of the single 
particle trajectory results. 



2.3. Test Particle Code and Parameter Choice 



We modify the t est p article code of iDalla fc Browning] 
(|2005l I2006L 120081) and iBrowning et all (|2010H to use the 
electromagnetic fi elds given above (from the solutions of 
ICraig et al1ll997f) . A Variable-Step Variable-Order Adam's 
method, where the step size is recalculated to properly 
resolve gyro-motion, is used to integrate the relativistic 
Lorentz equation. 



dp 
~df 



m 



E 



P 

7m 



x B 



(28) 



where p is the momentum of the particle, 7 is the Lorentz 
factor, q and m are the charge and rest mass and E and B 
are the analytic expressions for the electric and magnetic 
fields for each model. 

We use the expressions for the electric potential V, cal- 
culated in equations (|12|) and ((23)) to calculate the electric 
potential energy at each time step. With this we verify that 
the total energy 

W = e k + qV (29) 

is conserved where e& = (7— l)mc 2 . For each simulation we 
find that this is conserved up to 5 significant figures. Also, 
to check the code handles non-adiabatic motion in strong 
magnetic field gradients of a current sheet we reproduce the 
results of lSpeiseil (|1965l) . including the ejection time for the 
case with background field. 

We choose Lq, the normalising length scale, to be Lo — 
10 4 m for global simulations to keep integration times short. 
This size of simulation box can be considered as a local re- 
gion around the null at which the linear background field 
and flow in equation (J5J are good approximations. We use 
a larger value of Lq = 10 6 m for simulations where particles 
are initially distributed within the current sheet, as veloc- 
ities are typically much faster here. Note that a change in 
L also changes the value of 77 as given in equation (TT|). 
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Fig. 3: Angular distribution of protons from the null point 
for spine model at t = 1.6 x 10 6 T WjP , at which the en- 
ergy spectrum is steady state, for parameters A = 0.75, 
B s = 10, a = —0.1, r) = 10~ 6 . The initial distribution was 
Maxwellian at T = 86 eV in the upper right inflow region. 



All magnetic fields mentioned in Section [2] have dimen- 
sions of Bq = 0.01 T, typical for the solar corona. We set 
vq = VAe = 6.5 x 10 6 ms _1 (corresponding to a number 
density no = 1-126 x 10 15 m~ 3 ). All dimensionless times 
quoted are in terms of the gyro-period, T u — 2nm / '(qBo) . 

To examine single particle trajectories in both models 
we choose values of B s = 10, A = 0.75, r\ = 10~ 6 . This 77 
value is rather large, towards the highest possible anoma- 
lous resistivities (with Lq = 10 4 m), but is useful to ob- 
serve particles entering the current sheet. We vary all three 
of these parameters to produce scalings of energy gain and 
drift times, where we use values as low as those expected 
in purely collisional plasma i.e. rj = 10~ 12 . 



3. Results 

3.1. Spine Global Trajectories 

Initially, we place a distribution of 5 000 protons with 
Maxwellian velocities of temperature T = 10 6 K (86 
eV) in the spine model fields. The protons have positions 
from a uniform random distribution at a global distance 
x 2 + y 2 + z 2 = 1 from the null point. We only discuss here 
protons that start in the upper right inflow region of longi- 
tude < <t> < 180° and latitude < (3 < 90° (here (f> = is 
the x-axis and /3 = is the fan plane). 

Figure [3] shows the final spatial distribution and ener- 
gies of the particles at time t = 1.6 x 10 6 T UiP ~ 10 s, at 
which the energy spectrum in Figure [B] becomes steady- 
state. Those protons starting in the lower left inflow region 
have final distributions as in Figure [3] after reflections in 
both = and (3 = apart from statistical differences. The 
parameters used here are B s — 10, r\ = 10~ 6 , a = —0.1. 
This value of a limits the bulk flow exhaust speed to be 
sub-Alfvenic but it increases the electric drift speed in the 
external region (see equation (|19D ) due to weaker magnetic 
field on r — 1. This gives reasonable simulation times, but 
there are still some particles in the upper right inflow quad- 
rant at the end of the simulation. 



6 



A. Stanier et al.: Solar Particle Acceleration at Reconnecting 3D Null Points 




Fig. 4: Typical proton global trajectory from population 'A' 
for parameters A = 0.75, B s = 10, a = —0.1, r\ = 10~ 6 . 
The particle is taken from the many particle simulation 
having initial position x = (—0.52, 0.80, 0.29) and velocity 
v = (-0.0044, 0.0013, -0.0088)?M e . The magnetic field lines 
(thin dashed) are a projection of the field from the plane of 
the trajectory <j> ss 120°. Inset shows the 3D trajectory of 
the proton as it crosses the spine-axis, the solid line in the 
centre is the line B z (x, y, 0.95) = 0. 



There are two main populations of accelerated parti- 
cles. The population labelled 'A' in Figure |3] is close to the 
fan plane, \f}\ < 10°, with energy e& > 1 keV, and with 
longitude —90° < <j> < 90° comprising of about 8% of the 
total proton number. The maximum particle energy of this 
population is about 15 keV. Note that the current in the 
spine axis is aligned with <fi — through the centre of this 
population. There are also some high energy protons scat- 
tered at large positive latitudes for <j> < 0, and at large 
negative latitudes for (j> > 0. To look more closely at what 
is happening here we will choose a typical proton from this 
population and follow its trajectory below. 

For those particles that have crossed the fan plane, 
P = 0, into the lower right outflow quadrant, the spatial 
an d energy distributio n look s similar to the ideal spine case 
m iDalla k Br owning (2006). The accelerated population 
which has e k > 1 keV and (3 < -85° is labelled 'B'. This 
population is about 6% of the total protons in the simula- 
tion and the maximum kinetic energy in this population is 
£k,max ~ 12 keV. The angular distribution differs slightly 
with the ideal case in that there are few particles found 
between the latitudes —70° < /? < —85°; particles appear 
to be closer to the negative spine axis in the resistive case. 




-0.02 -0.01 0.00 0.01 0.02 

x 

Fig. 5: Typical proton trajectory from population 'B' 
in the many particle simulation, with initial posi- 
tion x = (—0.54, 0.78, 0.31)L and velocity v = 
(-0.004, -0.006, 0.002)v Ae . The dashed lines show the pro- 
jection of the magnetic field from the plane of motion, 
<f> 110°, onto the y-z plane. Inset shows the motion in 
the x-y plane close to the spine axis. The purple arrows 
show the direction and relative magnitude of the gradient 
drift velocity and the dash-dotted lines show contours of the 
electric potential, with the intersecting tick mark indicating 
lower potential to the right. 



A typical proton trajectory from population 'A' is 
shown in Figure|H The proton which starts at (xq, yo, zq) = 
(—0.52, 0.80, 0.29) in the upper right hand inflow quadrant 
initially moves away from the null but mirror bounces and 
travels back towards the spine along the fan plane. The elec- 
tric drift speed increases towards the spine causing the pro- 
ton to enter the resistive region, which has radius r n 0.01, 
about the spine axis. It enters at (x, y, z) ps (—0.01, 0, 0.95) 
after t = 3 x 10 5 T WjP « 2 s (inset). At this point the 
proton becomes unmagnetised as the gyro-radius becomes 
comparable to the length-scale of magnetic field gradient 
p/L^B > 1 (we typically find that gyro- motion starts to 
break down when p/Lys ^ 10~ 2 ). 

The proton is then directly accelerated in the x-direction 
parallel to the current at the spine with du/dt s» qEo/m as 
it crosses r = 0. For small displacements in the y-direction 
a strong Lorentz force due to the B z field returns it to y = 
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line. These oscillations are Speiser-like (|Speiserlll965l) with 
frequency approximately ui oc t 1 ! 2 . 

The Speiser-like motion finishes and the first gyrations 
start (not shown) when the proton reaches r ~ br^ at which 
p/LyB < 1- However, the energy gain of efe ps If keV is lo- 
calised to within x ~ 2r v , during which the trajectory does 
not deviate much from the x-direction (note the y-axis scale 
in the inset of Figure 0]). In effect, the proton has left the lo- 
calised current sheet while unmagnetised but before it can 
be ejected by the background field components, in contrast 
to 2D current sheet configuratio ns with weak guide field (eg. 
lSpeiser|[T965t lLitvinenkolll996ft . Figure H may give the im- 
pression that the particle is being ejected, however, this is 
just the centre of the Speiser-like oscillations following the 
B z (x,y,z — const.) = line (which here is not straight as 
in the usual 2D configurations). This behaviour is evident 
considering the F = qv x B force for the unmagnetised 
proton if B z is the dominant component of the magnetic 
field. 

After the proton becomes re-magnetised at r w 5r r) 
it has weak electric drift, ve <C v u . It follows the field- 
lines closely and mirror bounces travelling back towards 
the spine: there the proton is taken up to high latitude be- 
fore it bounces again. This mirror bouncing is the reason 
for the 'scattered' accelerated protons in Figure 02 some of 
which are at large latitudes. 

A typical particle trajectory chosen from population 'B' 
is shown in Figure^ The proton starts at (—0.54, 0.78, 0.31) 
and drifts towards the spine but bounces and crosses the 
fan plane instead. It exits the simulation box down the 
base of the spine axis, reaching an energy efe = 6.72 keV 
as it crosses z = —5. As there is no electric field in the 
z-direction, the energy gain must occur due to motion in 
the x-y plane, which is also shown in Figure [SI The proton 
enters the region close to the spine axis parallel to a contour 
of the electric potential, but then drifts across the contour 
due to strong gradient drift. While the proton gains energy, 
the gradient drift is larger than the electric drift by a factor 
of 2 with the latter directed inwards towards the current 
sheet. The proton is stopped as it reaches z = —5Lq which 
we do consistently throughout these simulations. At the 
time of stopping it is actually losing energy as it re-crosses 
the same electric potential contours. However, some other 
protons from the many-particle simulation in Figure|3]reach 
the current sheet at low latitudes, gaining higher energy. 

The energy spectrum for the spine simulation is shown 
in Figure El If protons cross the R = 5Lq spherical bound- 
ary we use the energy at the instant of crossing (if this 
is not done some protons reach order ~ 10 2 Lq which be- 
comes unrealistic as the background field increases with- 
out bound away from the null, also causing the time-step 
to decrease and simulation time to increase). The initial 
Maxwellian spectrum hardens to a broken power law with 
maximum energy of about €k ~ 15 keV. This maximum 
energy can be understood as the difference in potential en- 
ergy across the spine current sheet, ~ qEx acc where 
E ps Eq ss nB s /r n [vAeBo] and x acc ss 2r v [Lq] is the accel- 
eration distance (from —r v < x < r n ), as the electric field 
drops off quickly for \x\ > r v . For the parameters used, this 
gives e^. ps 13 keV. This approximate expression has no de- 
pendence upon the parameter a, so the limiting of a < B s 
should not have a large effect on this result. 
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Fig. 6: Energy spectrum from the many particle simulation 
for protons in the spine model, with parameters A = 0.75, 
B s = 10, a = —0.1, i] = 10 -6 . For particles leaving the 
R = 5 sphere, the energy at the time of crossing is used. 



3.2. Fan Global Trajectories 

The many particle simulation for the fan model is shown 
in Figure [7] for the optimised solution B s — a — 10, 
with T] = 10~ 6 , A = 0.75. The initial distribution has 
thermal energy e& = 86 eV with uniform random posi- 
tion in the upper inflow quadrant —90° < <f> < 90° and 
< f3 < 90°. The final angular distribution is taken from 
when the proton distribution reaches a steady state in en- 
ergy at t = 4000 T w p m 0.025 s. This is more than two 
orders of magnitude faster than the spine model for simi- 
lar parameters (even after the spine drift was increased by 
limiting a < B s ) as the external electric drift (equation |26| 
scales more favourably with the resistivity. Protons that 
cross the R = 5Lq spherical boundary from the null point 
before this time are stopped and the energy and angular 
position at time of crossing is used. The t = angular dis- 
tribution has some structure in terms of final energy gain 
as the initial random thermal velocities are dominated by 
the strong electric drift. 

Within this structure there is some asymmetry in (f>. 
Indeed, we do not expect symmetry between particles drift- 
ing clockwi se and anti-clockw i se ab out the null as in the 
ideal case (jDalla fc Browning 120061 ) now that there is a 
current in the y direction. Those protons with e/. ~ 10 7 
eV at <f) fa —20° (the yellow vertic al b and to the left of 
the green vertical band in Figure [ ffia)| do not enter the 
current sheet, but gain high energy, as they are unmag- 
netised slightly, p/L^B ~ 10~ 3 , due to very fast electric 
drifts close to the sheet. Here the first adiabatic invariant, 
the constancy of fi, is also violated. 

Typically, the high energy protons of Figure [7] start ei- 
ther close to the x-axis at low to mid latitudes (about 7% 
of the total number at latitude f3 > 1° with final energy 
ekjin ^ 10 MeV), or they start at very low latitude close 
to the fan plane (< 1% of total at /3 < 1° and Ckjin ^ 10 
MeV). At t = 4000 T WiP these energetic protons are found 
at (3 pa either side of 4> = 90°; the y-axis in the direction 
of the fan current. 



8 



A. Stanier et al.: Solar Particle Acceleration at Reconnecting 3D Null Points 



50 



■5 




-100 100 

longitude (<p) 



(a) 



S3. 

^ o 



-50 




-100 100 

longitude (0) 

(b) 

Fig. 7: a) Angular distributions of protons in fan model at 
t = 0, with initial temperature T = 86 eV. The x-axis is 
= and the fan plane is /3 = 0. Protons are coloured by 
the final energy at t = 4 000. Parameters used are A = 0.75, 
B s = a = 10, n = 10~ 6 . b) Angular distributions at time 
t = 4 000T U)P when the energy spectrum has reached steady 
state. 



At i = 4 000 T^p there are a small number of high en- 
ergy protons scattered at high latitudes (about 0.1% with 
efc > 10 MeV). These enter the current sheet temporarily 
at negative longitude far from the null point, but exit again 
without any Speiser-like motion. They become slightly un- 
magnetised, with maximum p/iys ~ 10~ 2 , following com- 
plicated trajectories. As they are not typical they are not 
investigated further in the external region, but the be- 
haviour within the current sheet is discussed below (shown 
in Figure IT21. 

Figure |S] shows the trajectory of two typical protons 
taken from the simulation. Proton T' starts at (xn, y®, zq) — 
(0.86, 0.41, 0.30) and drifts around the null point due to the 
strong azimuthal electric drift. Although it drifts down to- 
wards the current sheet, it reaches a minimum height of 
z ~ 15 z v before it flows into the outflow quadrant, not 
entering the sheet. The main velocity contribution is elec- 
tric drift as it moves around the null point, but v\\ becomes 
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(b) 

Fig. 8: Two typical proton trajectories from the many par- 
ticle simulations in the fan model. Proton T' is represented 
by a thin line with initial position (0.86, 0.41, 0.30), and pro- 
ton '2' by a thick line with initial position (0.8, 0.003, 0.6). 
The parameters are A = 0.75, B s = a = 10, rj = 10 -6 . The 
solid lines are representative magnetic field lines (seeded 
from the top of the spine axis and projected into the 2D 
planes) and the arrows show the direction and relative mag- 
nitude of the electric drift velocity, a) In the x-y plane, 
where the electric drift arrows are from the edge of the cur- 
rent sheet z = z v . b) In the x-z plane close to the current 
sheet, where the electric drift arrows are plotted on y = 0. 
The initial positions are not shown in this plane. 



dominant as the particle exits the simulation box paral- 
lel to the negative x-axis. The first adiabatic invariant is 
not violated, /i = const, and the maximum p/L\js ~ 10~ 4 
at closest point of approach to the sheet. The proton is 
strongly magnetised throughout. Despite not reaching the 
current sheet the energy gain is still considerable, reaching 
0.5 MeV as it crosses the R — h sphere. 
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Fig. 9: Energy spectrum for the many particle fan simula- 
tion for protons with parameters A = 0.75, B s = a = 10, 
77 = 10 -6 . For particles leaving the R = 5 sphere, the energy 
at the time of crossing is used. 



Particle 2 starts at (0.8, 0.003, 0.6)L . The azimuthal 
electric drift is weak close to the x-axis and the proton drifts 
down to the fan current sheet. It enters the sheet at (x, y) — 
(—0.02,0.15) and becomes unmagnetised: p/Lvs > 1 and 
fi is not conserved. We observe Speiser-like oscillations as 
the proton is accelerated in the y-direction. At t = 0.846 ms 
after entering the sheet, it passes out of the simulation box 
at R = 5. Here, the particle is still within the sheet with 
vn — 0.36c and = 67 MeV. Using this time period in the 
direct acceleration formula, y = qEot 2 /2m, with the electric 
field on z = 0, E = -nBs/f] 1 / 2 [v Ae B ] from ([22]). gives 
j/«5. Thus the proton is directly accelerated in the current 
sheet for the entire length of the simulation box. However, 
this motion is not Speiser-like throughout as p/Lvs < 10~ 2 
when the proton reaches y = 1.5Lrj- The proton reaches a 
global distance in the y-direction and becomes magnetised 
by the background B y component of the magnetic field, 
which acts as a kind of guide field. When the simulation 
is run without stopping the particle at R = 5, the proton 
is not ejected from the current sheet throughout the whole 
simulation time t = 4 000 T u p . 

This particle enters the current sheet at a distance 
R ~ 0.15 from the null point; however, this distance is 
not typical for the many particle simulation in Figure [7] 
In the simulation 9.3% of the total particles reach the cur- 
rent sheet, after a mean time of about 800T UiP . The average 
distance from the null point of particles entering the sheet 
is R rs 2.2; some remain magnetised by the background 
magnetic field inside the sheet. 

The energy spectrum for the fan simulation is shown in 
Figure [H Almost all the protons are accelerated into a non- 
thermal power law distribution f(E) oc E~ J , with slope 
7 sa 1.5. For most particles, this efficient acceleration is 
due to the fast electric drift speed in the fan model being 
much larger than the initial thermal velocity. The spectrum 
appears to have reached a steady state by t = 4 0007^; 
however, it also depends upon the position at which pro- 
tons are stopped as they leave the simulation box. As a 
test we repeat the simulation but stopping the protons at 



a spherical surface of radius R = 10 from the null point, 
instead of R = 5 that has been used consistently through- 
out these simulations. Now the 'fiat tail' at 10 75-8 eV in 
Figure O becomes a 'bump on tail' centred at 10 8 eV (not 
shown) disconnected from the main distribution. Here, the 
power law part remains mostly unchanged. The population 
of protons that is trapped in the sheet as it crosses R = 5 
due to the strong 'guide field' remains trapped at R = 10 
where B y (y) has doubled in strength. 

3.3, Fan Current Sheet Trajectories 

The simulations considered thus far concern proton tra- 
jectories starting from the external region, at a distance 
R = 1 from the null point. However, most of the protons 
entering the current sheet do so far from the null point. In 
the following, protons are initially distributed within the 
fan current sheet close to the null, to study the transition 
from non-adiabatic to adiabatic motion. 

Firstly, we place particles within the sheet so that they 
are initially unmagnetised by the B y {y) component of the 
background field. They are magnetised only by the strong 
B x (x, z). The protons are uniformly distributed in the area 
I a: I < 1; y = 0; \z\ < z v with initial thermal energy 
T = 86 eV. Figure [TD] shows the position of 2 000 protons at 
t = 2 500T WiP (a), and t = 17500T^ iP (b), during this sim- 
ulation for the parameters A = 0.75, B s = a = 5, r\ = 10 -8 . 
We increase the dimensional box length to Lq = 10 6 m as 
velocities in the current sheet are typically fast, giving rea- 
sonable integration times. This makes our result s more com- 
parab le to the approximate analytic solutions of lLitvinenkol 
(2006J). Note that rj decreases due to the increase in Lq in 
equation (JXJ) . We again artificially stop the particles as they 
cross the R — 5 spherical surface. 

At t = 2 500 T up most of the protons are strongly mag- 
netised by the B x (x,z) magnetic field. Inside the current 
sheet, \z\ < z n we can use equation (JT3J) to get approximate 
expressions for the electric and magnetic fields, 



E y y « r)B s /z v y [v Ae B ], 



B 



( Xax 



B s z Xay 



-Xaz 



[B 



(30) 



(31) 



Zff, on the 



E z is small except at global distance in y (see below). 

For a proton starting at x — 0, y — 0, z 
edge of the current sheet, the background components of 
the magnetic field are negligible. The proton drifts towards 
the vertical centre of the sheet ve z ~ [vAe]- It 

becomes unmagnetised at the fan plane, z ~ 0, close to the 
null point and is directly accelerated in the y-direction. We 
compare th i s traj ectory to the analytical WKB solutions of 
iLitvinenkol ((2006). The ejection time for a non-relativistic 
proton that is unmagnetised close the null point, x ~ 0, z as 
0, in the fan current sheet in our parameters is 



L ejec 



m 2 BeL l 



q 2 z v B X 2 a 2 Ey 



1/3 



(32) 



(jLitvinenk o 2006) provided that the proton remains within 
the non-adiabatic region and the displacement magnetic 
field gradient is much stronger than the gradient from the 
background component, B s /z v ^> Ace. The second assump- 
tion is valid for our simulation; however, we do not observe 
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(a) t = 2 500 




(b) t = 17 500 

Fig. 10: Proton positions after initial distribution within the 
fan current sheet, such that |xo| < 1) \yo\ = 0, \zq\ < z v = 
y/2fj. Parameters used are A = 0.75, B s = a = 5, rj = 10~ 8 , 
Lq = 10 6 m. The dashed lines are representative magnetic 
field-lines inside the current sheet (note the difference in 
scale of the z-axis). The solid black line is the line (x±, 0, z\) 
such that B x (x\, 0, Z%) = 0. The particles are stopped at 
R = 5 from the null point. 



proton energy gain limited by ejection in these simulations. 
To understand this, we consider the distance travelled in 
the y-direction during this time, 



Veje 



y{te 



qE y t 2 



'ejec 



2m 



(33) 



which we compare with size of the non-adiabatic region 
from the null in this direction. The particle begins to be 
re-magnetised by the background field at a global distance 
y* such that 

v{y')/y' * Wflvto .) (34) 

where v(y) is a typical proton velocity and UB y (y) is the 
gyro- frequency of a particle gyrating around B y (y). We use 

1/2 

v(y) = (2qE y y/m) 1 from direct acceleration (if we use 
v(y) = Ey/By the value for y* differs by 2 1 / 3 ), assuming 
that there was no initial y- velocity and the particle entered 



the sh eet at y rj 0. We recover the result of iLitvinenkol 
(|2006f h that in dimensional form 



SmE, 



v 



q(B Xa) 2 L 



The ratio of these two distances is 
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ejec 
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(35) 



(36) 



where we have ignored factors of order unity. The ratio of 
the two timescales is the square root of this. There is little 
gyro-turning for protons starting close to the null point 
as this ratio is necessarily small for the fan current sheet. 
The proton is magnetised by the B y (y) "guide field" and 
trapped in the sheet, the energy gain is only bounded by 
the length of the sheet. 

Figure QJJ] also shows the more general case of protons 
starting at y = 0, \z\ < z n and at a global distance in x. 
These protons drift vertically until they reach the diagonal 
line where B x {x,z) = 0, at which they become unmag- 
netised and accelerated. We do appear to see some gyro- 
turning for protons starting at \x\ rs 1. This is probably 
due to the strong component of the Lorentz force, v y B z , 
that acts to turn the trajectory to the x-direction. For par- 
ticles starting at x — the B z magnetic field switches sign 
during the z-oscillations, but at x = 1 the proton is un- 
magnetized below the fan plane and B z stays positive. The 
proton is turned in the positive x-direction but is quickly 
magnetised by the guide field when it reaches a distance of 
about y* (see equation 133)) . In Figure [TU] it can be seen that 
particles are accelerated radially outwards from the null. 
They continue to gain energy as they become magnetised 
about the background field P, on a field-line with a parallel 
component of the electric field En = E ■ P/\P\. 

We artificially stop the protons at R = 5 Lq from the 
null point. At t = 50 000 T u , p all of the protons in the simu- 
lation have reached this distance without being ejected and 
we fit the energies of the particles by the expression 



«5 9 ?7 1/2 B 3 J 2 (1 



X 2 ) 1 ' 2 sin^ [vAeB L ], (37) 



for the optimised solution a — B s , where (f> is the azimuthal 
angle (4> = 90° is parallel to the current). Figure ITTI shows 
the energies of 5 000 protons in three simulations with iden- 
tical setup to Figure [TUJ but with different values of r\ and 
B 8 . This expression (thin line) fits the energies of simulated 
particles (circles) as they cross R = 5 Lq very well. 

In Figure [T^] we place 5 000 protons in the fan current 
sheet with initial position in \z\ < z r) , —1 < x, y < 1 so that 
they are initially magnetised by the "guide field" B y (y). 
This is the more general case, as protons reaching the cur- 
rent sheet from the external region will not typically do so 
at y sa 0. The protons that do not start close to y = 
are directly accelerated without the initial drift phase. By 
t = 19 000 T^p all of the protons have left the simulation 
box; either through the R = 5 boundary, or through the 
edge of the current sheet \z\ — z v . The particles that cross 
\z\ = z n in y > start close to the edge and leave due 
to initial thermal velocity. However, those starting with 
y < are ejected from well within the current sheet. These 
protons (19.7% of total number) are circled in Figure [P2l 
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Fig. 11: Energy distribution of particles as they cross the 
R = 5 boundary, where initial position is within |a;o| < 1, 
\yo\ = 0, \z \ < z v = y/2f]- Here, L = 10 6 m and the results 
for different values of r\ and B s are plotted. The solid points 
are protons from the three simulations and the thin lines 
show the suicj) relationship in equation (|37[) . 



Typically they remain magnetised, with p/Lys m the range 
10 -4 — 10~ 2 . They are not ejected due to gyro-turning in 
the sense of Spciscr (1965) as this requires non-adiabatic 
motion. The trajectories in the y-z plane seem to follow 
the magnetic field lines closely, although they have strong 
electric drift from the E z component of the electric field. 
Within the current sheet, \z\ < z n the truncated power se- 
ries in equation (|13[) gives the z-component of the electric 
field from equation ([^1 as 



B s a 



(1 - X 2 )yz, 



(38) 



which is stronger than the current electric field (|3T)|) for 
global y and 0. However, it only contributes to strong 
electric drift (not shown) in negative x-direction for protons 
in the upper half of the sheet < z < z n , and in the 
positive x-direction for —z n < z < 0. This electric field 
also contributes to ve v but this is dominated by the direct 
electric field acceleration. 

The protons that are not ejected from the current sheet 
have an approximate sinusoidal dependence in kinetic en- 
ergy gain, given by equation (|37p (there is a thicker spread 
of points about the predicted lines than in Figure [TT] due 
to differences in initial potential energy). 

3.4. Scalings 

In the fan global simulation of Figure [7] we chose the opti- 
mised parameters r\ = 10~ 6 , B s = a = 10 and A = 0.75. 
With consideration to the large variation in both scale and 
behaviour in a given distribution of flares, it is interesting 
to see how the results of this simulation scale when the 
simulation parameters are varied. 

The current sheet within the fan model is the most effec- 
tive way to accelerate the particles. Thus, it is interesting 
to study the effect of varying parameters on the fraction 
of protons that enter the current sheet from the external 
region, and the average time taken to drift there from an 
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Fig. 12: Proton positions after initial distribution within the 
fan current sheet such that \xo\, \yo\ < 1 an d \zq\ < z v — 
\[ T £f\. Parameters used are A = 0.75, B s = a = 5, r\ = 10 -8 , 
Lq = 10 6 m. Particles circled in black are those that start in 
y < and cross z = z v « 9.6 x 10~ 5 before t = 19 000 T UtP . 



initial position on the R = 1 sphere. These scalings are 
shown in Figure 1131 They are from simulations of 5 000 
protons at T = 1 MK starting at the upper inflow region 
at R = 1. We define the current sheet as z — z v — for 
the fan model, although we note that not all of the protons 
reaching this height become non-adiabatic. 

The average time taken for the particles to reach the 
current sheet gives a measure of the external electric drift 
speed. The approximate drift scaling of equation (|2l)|). ve ~ 

Bs^r) 1 / 4 , is in reasonable agreement with these drift times. 

The fraction of particles reaching the sheet increases 
typically with increasing r\ and decreasing B s . Note that 
the size of the current sheet which we use to produce the 

scalings has the dependence z v ~ t^I 2 B s ^ 2 (1 — A 2 ) -1 / 2 , 
although this does not fully explain the result as the trends 
are not simple power laws. Figure PHI shows how the energy 
spectra vary with these parameters. As might be expected 
the spectra shift to the right for an increase in both B s and 
r\. Both the convective electric field (and so external electric 
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Fig. 13: Percentage of total particles (+) reaching current 
sheet at z = from the external region {R = 1 in the inflow 
quadrant), and the mean time taken (*), for different values 
of rj, B s and A. Each data point is from a many particle 
simulation with initial Maxwellian distribution (T = 86 
eV) of 5 000 protons. The set-up is the same as that in 
Figure [7] a) For different values of r] with fixed B s = a = 5, 
A = 0.75. The solid line is a least squares fit to the points, b) 
For different values of B s (with B s — a) for fixed rj = 10~ 8 , 
A = 0.75. The solid line is a least squares fit. c) For different 
A with fixed B s = a = 5, rj = 10~ 8 . No curve was fit. 
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Fig. 14: Scalings of steady-state energy spectra for global 
fan simulation, a) For different values of rj with fixed B s = 
a = 5 and A = 0.75. The time taken to reach steady state 
was t = 8 x 10 4 , t = 2 x 10 4 and t = 6.4 x 10 3 T UtP for 
X] = 10~ 10 , r\ = 10~ 8 and r\ = 10~ 6 respectively, b) For 
different values of B s — a with fixed 77 = 10 -8 and A = 0.75. 
Steady-state was reached at t — 2 x 10 5 , t = 2 x 10 4 and t = 
5 x 10 3 T u , p for B s = 1, B s = 5 and B s = 30 respectively, 
c) For different values of A with fixed B s = a = 5 and 
r] = 10~ 8 . Time to steady state is t = 10 5 , t = 2 x 10 4 
and t = 8 x 10 3 T u p for A = 0.9, A = 0.75 and A = 0.3 
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drift) and the direct electric field within the sheet increase 
with larger values of these parameters. 

Up until now we have used A = 0.75 as a constant in all 
of the simulations. This parameter has been typically left 
as const ant in the calculation of MHD energy dis sipation 
scalings (|Craig et all 119971; ICraig fe Watson! |2000D for the 
fan model as it can only be varied within an order one 
range. However, it has a large effect on the efficien cy of the 
fan model for particle acceleration (see Figure ll^c)| and 
Figure fT^c)| . Varying A within < A < 1 has a comparable 
effect on the fraction of particles reaching the current sheet 
as varying rj by six orders of magnitude. Also, as shown in 
FigureHH decreasing A shifts the energy spectrum to higher 
energy and decreases the time taken to reach steady-state. 
These effects can be explained somewhat by an increase in 
external drift speed, although varying A also has an effect 
on other quantities such as the current sheet height. 



4. Summary and Discussion 

We investigat ed test part i cle m o tion in the elect r omag - 
neti c fields of ICraig et al.l (fl995h : ICraig fe Fabling! ([19961) 



and ICraig et aL ~ (|l997l) . that are solutions to the steady 
state, incompressible and resistive MHD equations at a 
3D null point. The study was carried out by mo dify- 
ing the code of iDalla fe Browning! (|2005l 120061 ^008) and 
iBrowning et all (|2010[ ). We considered initially Maxwellian 
(T = 86 eV) distributions of protons starting at a global 
distance R = L$ from the null point in resistive spine recon- 
nection, where the electric current is within a thin cylinder 
about the spine axis, and resistive fan reconnection, with a 
current sheet in the fan plane. When the energy spectrum 
from the simulations reached steady-state we find the final 
angular position of the particles from the null and their 
energy distribution. We identified different populations of 
accelerated particles and, to understand the acceleration 
mechanism, examined a typical single particle trajectory in 
each case. For the fan model we ran additional simulations 
with the particles initially distributed within the fan cur- 
rent sheet, to study the effect of the null point on directly 
accelerated particles. We consider two cases, where particles 
arc firstly unmagnetised and secondly magnetised in their 
initial position by a background "guide field". Finally, we 
show how the external drift times and energy spectra for 
the fan model scale when treating A, B s and rj as free pa- 
rameters (for the optimised so lution a = B s , correspo nding 
to the thinnest current sheet (jCraig fe Watsonll2000f) l 

We found that the spine model, wh ich gave promising 
accel eration results in the ideal case (jDalla fc Browning! 
2008), is much less effective when resistive effects are in- 
cluded (at least for this specific resistive model). The elec- 
tric drift in the external region is weak, scaling with resistiv- 
ity as ve{t 3> r n ) ~ 77, giving very long drift times for pro- 
tons to reach the spine axis. We find that there are two pop- 
ulations of accelerated particles. One of these escapes the 
simulation box down the base of the sp ine axis, similar to 
the p roton jet found in the ideal model (jDalla fc Browning] 
2006), and the other is close to the fan plane, where par- 
ticles have crossed the current sheet in the spine axis. The 
energy gain for particles that reach the current sheet is 
low. The main limiting factor is the small electric potential 
energy difference across the current sheet, due to localisa- 
tion of the reconnection electric field to a small cylinder 
about the spine axis. The apparent contrast with the ideal 



results arises from parameter choice. In the ideal regime, 
the magnitude of the electric field for the spine and fan 
models was set to equal strength in the external region, 
at a global distance from the null. The electric field falls 
steeply as 1/r in the external region of the spine model 
which gives very strong acceleration close to the spine in 
the ideal case. However, when the pressure constraints in 
the resistive model are taken into account, namely the lim- 
iting of the displacement field on the edge of the current 
sheet to avoid unphysical magnetic pressures (jCraig et al.l 
Il997f ), this 1/r dependence gives very weak electric field, 
and slow drift, in the external region. 

We found much higher proton energies in the resis- 
tive fan model for similar parameters. For 77 = 10 -6 , 
B s = a = 10 and A = 0.75 we find the energy spectrum 
from a distribution of protons starting with thermal energy 
at R = Lq from the null becomes power law at steady state, 
with a spectral index of about —1.5 and maximum particle 
energy of the order 0.1 GeV. The electric drift in the exter- 
nal region is much quicker than the spine model, ve ~ rj 1 ^ 4 . 
It accelerates all of the particles in the simulation as it 
is faster than the initial random velocities associated with 
thermal motion at T = 86 eV. We find that the popula- 
tion with the highest energy gain corresponds to protons 
that have entered the fan current sheet. The energy gain 
for these protons is not limited by ejection due to unstable 
motion as they are re-magnetised within the current sheet 
by a "guide field" . The upper bound in energy gain is only 
limited by the electric potential energy, determined by the 
length of the current sheet. However, we find that a number 
of protons that enter the current sheet upstream of the null 
point can be ejected while remaining magnetised. This is 
due to the geometry of the background field lines, namely 
that they diverge at the null. We wi ll study this effect in th e 
futu re when we consi der electrons. IBrowning et al.l (|2010t ). 
and iGuo et all (|2010h show that electrons remain magne- 
tised at a closer distance to the null point, which may give 
a difference between the number of electrons and protons 
ejected in this manner. 

We find that the parameter A, which gives the degree of 
shear between the magnetic and velocity fields (such that 
< A < 1) has a large effect on the final energy spectrum of 
protons in the fan model. In the limit of A = the magnetic 
field in the fan model is annihilated. As we expect magnetic 
field to still exist in the reconnection site after a topological 
change it would be more likely that Awl. 

In these simulations we have neglected the electromag- 
netic effects of the non-thermal particles onto the back- 
ground fields. This is typical of the test particle approach, 
where it is assumed that the number of particles in the 
current sheet is a small fraction of the total number. For a 
large range of parameters in the fan model (see Figure 1131) 
this fraction is typically less than 5% of the total num- 
ber of particles starting in the inflow region. To estimate 
the strength of the magnetic fields from these particles, 
and the polarisation electric field from any charge seper- 
ation, it is necessary to also consider electrons which we 
will do in future work. A fully self-consistent approach, eg. 
using Particle In Cell simulations, is computationally ex- 
pensive at present, particularly in fully 3D geometries due 
to the large dynamic range of spatial scales. We have also 
neglected compres sibility, a simplifi cation used to get the 
analytic solutions (jCraig et allll997D . It would be interest- 
ing in future to include both time-dependence and resistiv- 
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ity, by using electromagnetic fields from MHD simulations, 
to see how particles behav e in so called spine-fan reconnec- 
tion (|Pontin et alj|2007afbh . 
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